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Abstract. The version 2.0 of the program SecDec is described, which can be used for the 
extraction of poles within dimensional regularisation from multi-loop integrals as well as phase 
space integrals. The numerical evaluation of the resulting finite functions is also done by the 
program in an automated way, with no restriction on the kinematics in the case of loop integrals. 



1. Introduction 

Nowadays we are in the long awaited situation of being confronted with a wealth of high energy 
collider physics data, enabling us to explore physics at the TeV scale. However, for the analysis 
and interpretation of these data, precise theory predictions are mandatory. In most cases, 
this means that calculations beyond the leading order in perturbation theory are necessary. 
Such calculations either involve integrations over loop momenta for the virtual corrections, or 
integrations over phase spaces for the real radiation corrections. In both cases, multi-dimensional 
parameter integrals need to be evaluated, which can contain ultraviolet singularities and, in the 
presence of massless particles, contain soft and/or collinear singularities. These singularities, if 
regulated by dimensional regularisation, appear as poles in 1/e, but factorising the poles from 
complicated multi-loop or multi-parameter integrals is a highly non-trivial task. The program 
SecDec, presented in [HE]) performs this task in an automated way, based on the algorithm 
of sector decomposition [31 HI [5]. Other implementations of sector decomposition in public 
programs can be found in [6l[71[8]. However, the latter programs, including SecDec 1.0 [I], are 
not designed to cope with kinematics where physical thresholds can be present in the integration 
region. The program SecDec 2.0 [2j is able to achieve this task, by an automated deformation 
of the integration contour into the complex plane [£l [10] \TT\ [12] . In this talk the emphasis is on 
the presentation of the new features of the program SecDec. 



2. The algorithm 

2. 1 . General framework 

The decomposition algorithm is described in detail in [13j . here we only describe the basic 
concepts. Consider an iV-dimensional parameter integral, for example 

In = [ dxi . . . [ dx N x^ 1 " 6 ,if ^ tt^tt , (1) 

Jo Jo [xi gi{x) + x 2 g2{x)\ 

1 Based on a talk given at the 14th International Workshop on Advanced Computing and Analysis Techniques 
in Physics Research (AC AT), Uxbridge, London, UK, September 2011. 



where 51,52 ar e polynomial functions which do not vanish for x%,X2 — > 0. However, we would 
like to factorize the integrand such that the term in square brackets is non- vanishing in the limit 
x\ — > 0, X2 — > 0. This can be achieved by a decomposition into sectors where x\ and X2 are 
ordered: we multiply equation (pQ) with 1 = Q{x\ — X2) + 0(^2 — x\) and substitute X2 — > x\ X2 
in the first sector and x\ — > X2 x\ in the second sector, to arrive at 

I N = [ dxi... f dx N x^ £ { ^ 3 - + x 2 l ~ e - — — ^— — —X , (2) 

Jo Jo I [01 0*0 + x 2 g2{x)\ [x igi {x)+g2{x)]j 

where (ji,g~i are functions of the transformed variables. As 51, 52 are non- vanishing at the origin, 
the terms in square brackets cannot lead to singularities anymore; the singularities are factored 
out into the monomials of type x~ l ~ € instead, and subtraction of the poles and expansion in e 
are straightforward in this form. In more complicated cases, the decomposition procedure may 
have to be iterated to achieve a full factorisation. A detailed description can be found e.g. in 

U S3]. 

2.2. Loop integrals 

A scalar Feynman integral G in D dimensions at L loops with N propagators, where 
the propagators can have arbitrary, not necessarily integer powers Uj, has the following 
representation in momentum space: 



g - /r[d 



N 

n^(W,M,m|; 



d% = ^ d D k , Pj({k}, {p}, mf ) = q) - ra) + id , (3) 

where the qj are linear combinations of external momenta pi and loop momenta ki. Introducing 
Feynman parameters leads to 

(_l)N„ °f N N u N v -{L+l)D/2 N 
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where U is a polynomial of degree L and J- is of degree L + 1 in the Feynman parameters. The 
Lorentz invariants which can be formed from the external momenta of the diagram, as well as 
propagator masses, are contained in T . As a simple example, consider the function T for a 
massless one-loop box diagram: 

T = — S12 X2X4 — S23 xix 3 — iS . (5) 

The functions U and J- can also be constructed directly from the topology of the corresponding 
Feynman graph I15j. and the implementation of this construction in SecDec version 2.0 is 
one of the new features of the program. 

For a diagram with massless propagators, none of the Feynman parameters occurs 
quadratically in the function T = J-q . If massive internal lines are present, J- gets an additional 

N 

term F{x) = Fo{x) + U(x) ^ XjUVj. U is a positive semi-definite function. A vanishing U 

i=i 

function is related to the UV sub divergences of the graph. In the region where all invariants 
formed from external momenta are negative, which we will call the Euclidean region in the 



following, T is also a positive semi-definite function. Its vanishing does not necessarily lead to 
an IR singularity. Only if some of the invariants are zero, for example if some of the external 
momenta are light-like, the vanishing of T may induce an IR divergence. Thus it depends on 
the kinematics and not only on the topology (like in the UV case) whether a zero of T leads to 
a divergence or not. The necessary (but not sufficient) conditions for a divergence are given by 
the Landau equations [16] : 

Xj(qj-m?) = Vj (6) 

W^ Xj (^)-™')=0. (7) 
1 j 

If all kinematic invariants formed by external momenta are of the same sign, the necessary 
condition T = for an IR divergence can only be fulfilled if some of the parameters Xi go to 
zero. These singularities can be regulated by dimensional regularisation and factored out of 
the function T using sector decomposition. The same holds for dimensionally regulated UV 
singularities contained in IA. However, after these singularities, appearing as poles in 1/e, have 
been extracted using sector decomposition, for non-Euclidean kinematics we are still faced with 
integrable singularities related to kinematic thresholds. How we deal with these singularities 
will be described in Section 14.41 

2.3. General 'parameter integrals 

The program can also deal with parameter integrals which are more general than the ones 
related to multi-loop integrals, for example phase space integrals involving massless particles, 
where the regions in phase space corresponding to soft and/or collinear configurations lead to 
singularities which can be extracted as poles in 1/e. The only restrictions are that the integration 
domain should be the unit hypercube, and the singularities should be only endpoint singularities, 
i.e. should be located at zero or one. We assume that the singularities are regulated by non- 
integer powers of the integration parameters, where the non-integer part is the e of dimensional 
regularisation or some other regulator. The general form of the integrals is 

„i „i m 

1= dx t ... dx N T\Pi(x, {a}P , (8) 
Jo Jo i=l 

where Pi(x, {a}) are polynomial functions of the parameters Xj, which can also contain 
some symbolic constants {a}. The user can leave the parameters {a} symbolic during the 
decomposition, specifying numerical values only for the numerical integration step. This way 
the decomposition and subtraction steps do not have to be redone if the values for the constants 
are changed. The V{ are powers of the form V{ = a% + fye (with ai such that the integral is 
convergent). Note that half integer powers are also possible. 

3. The SecDec program 

The program consists of two parts, an algebraic part and a numerical part. The algebraic part 
uses code written in Mathematica [T7] and does the decomposition into sectors, the subtraction 
of the singularities, the expansion in e and the generation of the files necessary for the numerical 
integration. In the numerical part, Fortran or C++ functions forming the coefficient of each term 
in the Laurent series in e are integrated using the Monte Carlo integration programs contained 
in the Cuba library [13 [19], or Bases [20]. The different subtasks are handled by perl scripts. 
The flowchart of the program is shown in Fig. [I] for the basic flow of input/output streams. 
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Figure 1. Flowchart showing the main steps the program performs to produce the result files. 
In each of the subdirectories loop or general, the file Template.m can be used to define the 
integrand. The produced files are written to a subdirectory created according to the settings 
given in param.input. By default, a subdirectory with the name of the graph or integrand is 
created to store the produced functions. This directory will contain subdirectories according to 
the pole structure of the integrand. 



The directories loop and general have the same global structure, only some of the individual 
files are specific to loops or to more general parametric functions. The directories contain a 
number of perl scripts steering the decomposition and the numerical integration. The scripts 
use perl modules contained in the subdirectory perlsrc. 

The Mathematica source files are located in the subdirectories src/deco (files used for 
the decomposition), src/subexp (files used for the pole subtraction and expansion in e) and 
src/util (miscellaneous useful functions). The documentation, created by robodoc^T\ is 
contained in the subdirectory doc. It contains an index to look up documentation of the source 
code in html format by loading masterindex.html into a browser. 

In order to use the program, the user only has to edit the following two files: 

• param. input: (text file) 

specification of paths, type of integrand, order in e, output format, parameters for numerical 
integration, further options 

• Template .m: (Mathematica syntax) 

— for loop integrals: specification of loop momenta and propagators, resp. of the topology; 
optionally numerator, non-standard propagator powers, space-time dimensions 

— for general functions: specification of integration variables, integrand, variables to be 
split 

The program comes with example input and template files in the subdiretories loop/demos 
respectively general/demos, described in detail in pQ. 

4. Installation and usage 

4-1- Installation 

The program can be downloaded from 
|http: //projects .h epforge . org/secdec, 

Unpacking the tar archive via tar xzvf SecDec. tar. gz will create a directory called SecDec with 
the subdirectories as described above. Then change to the SecDec directory and run ./install. 

Prerequisites are Mathematica, version 6 or above, perl (installed by default on most 
Unix/Linux systems), a Fortran compiler (e.g. gfortran, ifort), or a C++ compiler if the C++ 
option is used. 

4-2. Usage 

(i) Change to the subdirectory loop or general, depending on whether you would like to 
calculate a loop integral or a more general parameter integral. 

(ii) Copy the files param . input and Template. m to create your own parameter and template 
files myparamf ile, mytemplatef ile. 

(iii) Set the desired parameters in myparamf ile and specify the integrand in mytemplatef ile. 

(iv) Execute the command ./launch -p myparamfile -t mytemplatefile in the shell. 

If you omit the option -p myparamfile, the file param. input will be taken as default. Like- 
wise, if you omit the option -t mytemplatefile, the file Template. m will be taken as default. 
If your files myparamfile, mytemplatefile are in a different directory, say, myworkingdir, 
use the option -d myworkingdir, i.e. the full command then looks like ./launch -d my- 
workingdir -p myparamfile -t mytemplatefile, executed from the directory SecDec/loop or 
SecDec/general. 

(v) Collect the results. Depending on whether you have used a single machine or submitted 
the jobs to a cluster, the following actions will be performed: 



• If the calculations are done sequentially on a single machine, the results will be collected 
automatically (via results.pl called by launch). The output file will be displayed 
with your specified text editor. 

• If the jobs have been submitted to a cluster, when all jobs have finished, execute 
the command ./results. pi [-d myworkingdir -p myparamfile] . This will create the files 
containing the final results in the graph subdirectory specified in the input file. 

(vi) After the calculation and the collection of the results is completed, you can use the shell 
command ./launchclean [graph] to remove obsolete files. 

It should be mentioned that the code starts working first on the most complicated pole 
structure, which takes longest. This is because in case the jobs are sent to a cluster, it is 
advantageous to first submit the jobs which are expected to take longest. 

4 ■ 3. New features 

Version 2.0 of SecDec contains the following new features, some of which will be illustrated by 
examples in Section More details are given in [2] . 

• Multi-scale loop integrals can be evaluated without restricting the kinematics to the 
Euclidean region. 

• The possibility to loop over ranges of parameter values is automated, with the option of 
outputting results in a format suitable for plotting. 

• The user can define functions at a symbolic level and specify them only later after the 
integrand has been transformed into a set of finite functions for each order in e. 

• The regulator of the parameter integrals can be different from the dimensional regulator 
e. This is particularly useful to define e.g. measurement functions at a later stage of the 
calculation. 

• For scalar multi-loop integrals, the integrand can be constructed directly from the topology 
of the diagram. The user only has to provide the labels of the vertices connected by the 
propagators and the propagator masses, but does not have to provide the momentum flow. 

• The files for the numerical integration of multi-scale loop integrals are written in C++ 
rather than Fortran. For integrations in Euclidean space, both the Fortran and the C++ 
versions are supported. 

• Both the algebraic and the numerical part allow full parallelisation. 
4-4- Implementation of the contour deformation 

Unless the function T in equation ^ is of definite sign for all possible values of invariants and 
Feynman parameters, the integrand of a multi-loop integral will vanish within the integration 
region on a hypersurface given by the solutions of the Landau equations ([B]),©. However, we 
can avoid the poles on the real axis by a deformation of the integration contour into the complex 
plane. As long as the deformation is in accordance with the causal i5 prescription of the Feynman 
propagators, and no poles are crossed while changing the integration path, we can make use of 
Cauchy's theorem to choose an integration contour such that the integration is convergent. 
The i5 prescription for the Feynman propagators tells us that the contour deformation into the 
complex plane should be such that the imaginary part of T should always be negative. For real 
masses and Mandelstam invariants Sij, the following Ansatz [9j 110 1 ITT j 112] is convenient: 




x — i t(x) 




N-l 



(9) 



In terms of the new variables, we thus obtain 



Tim) = m - * a y, + °^ 2 ) > 




(10) 



such that T acquires a negative imaginary part of order A. 

The convergence of the numerical integration can be improved significantly by choosing an 
"optimal" value for A. Values of A which are too small lead to contours which are too close 
to the poles on the real axis and therefore lead to bad convergence. Too large values of A can 
modify the real part of the function to an unacceptable extent and could even change the sign 
of the imaginary part if the terms of order A 3 get larger than the terms linear in A. This would 
lead to a wrong result. Therefore we implemented a four-step procedure to optimize the value 
of A, consisting of 

• ratio check: To make sure that the terms of order A 3 in equation ()10p do not spoil the 
sign of the imaginary part, we evaluate the ratio of the terms linear and cubic in A for a 
quasi-randomly chosen set of sample points to determine a maximal allowed A = A max . 

• modulus check: The imaginary part is vital at the points where the real part of J- is 
vanishing. In these regions, the deformation should be large enough to avoid large numerical 
fluctuations due to a highly peaked integrand. Therefore we check the modulus of each 
subsector function T% at a number of sample points, and pick the fraction of the value of 
Xmax maximising the function with the minimal modulus, i.e. the value of lambda which 
keeps J- furthest from zero. 

• individual adjustments: If the values of are very different in magnitude, it can 
be convenient to have an individual parameter A(i, j) for each subsector function J-{ and 
each Feynman parameter Xj. 

• sign check: After the above adjustments to A have been made, the sign of Im^J 7 ) is again 
checked for a number of sample points. If the sign is ever positive, this value of A is 
disallowed. 

5. Examples and results 

In this example, we will demonstrate three of the new features of the SecDec 2.0 program: the 
construction of J-,IA directly from the topology of the graph, the evaluation of the graph in the 
physical region, and how results for a whole set of different numerical values for the invariants 
can be produced and plotted in an automated way. We will use the two-loop diagram shown 
in Fig. [2] as an example. Numerical results for this diagram have been produced in |22[ |23| . 
analytical ones in [241 12"5] . 




Figure 2. Two-loop vertex graph called P\2Q, containing a massive triangle loop. Bold lines 
are massive, dashed lines are massless. The vertices are labeled to match the construction of 
the integrand from the topology as explained in the text. 
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Figure 3. Comparison of analytic and numerical results for the diagram P\2Q- 



5.1. Topology-based construction of the integrand 

The template file templateP126 .m in the demos subdirectory contains the following lines: 
proplist = {{ms[l}, {3, 4}}, {ms[l}, {4, 5}}, {ms[l], {5, 3}}, {0, {1, 2}}, {0, {1, 4}}, {0, {2, 5}}}; 
onshell = {ssp[l] — s- 0,ssp[2] ->• 0, ssp[3] -> sp[l,2]}; 

where each entry in proplist corresponds to a propagator of the diagram; the first entry is the 
mass of the propagator, and the second entry contains the labels of the two vertices which the 
propagator connects. Note that if an external momentum p^ is flowing into the vertex, the 
vertex also must have the label k. For vertices containing only internal propagators the labeling 
is arbitrary. The on-shell conditions in the above example state that p\ = p\ = 0, p\ = s. 

5.2. Results in the physical region 

To run this example, from the loop directory, issue the command . /launch -d demos -p 
paramP126. input -t templateP126 .m. The timings for the finite part and a relative accuracy 
better than 1%, using Cuba-3.0 [19J, are about 100 seconds per point on an Intel(R) Core i7 
CPU at 2.67GHz, where the timings are very similar far from the s = 4m 2 threshold and close 
to threshold. 

5. 3. Producing data files for sets of numerical values 

To loop over a set of numerical values for the invariants s and m 2 once the C++ files are created, 
issue the command 

. /multinumerics .pi -d demos -p multiparamP126 . input. This will run the numerical 
integrations for the values of s and m 2 specified in the file demos/multiparamP126 . input. 
The files containing the results will be found in demos/21oop/P126, and the files p-2.gpdat, 
p-l.gpdat and pO.gpdat will contain the data files for each point, corresponding to the 
coefficients of e -2 , e _1 and e° respectively. These files can be used to plot the results against the 
analytic results using gnuplot. This will produce the files P126R0.ps, P126I0.ps which will 
look like Fig. [3j 

6. Conclusions 

We have presented the program SecDec 2.0, which can be used to factorise dimensionally 
regulated singularities and numerically calculate multi-loop integrals in an automated way. As a 
new feature of the program, it now can deal with physical kinematics, i.e. is not restricted 



to the Euclidean region anymore. A new construction of the integrand, based entirely on 
topological rules, is also included, and the new features are demonstrated for the examples 
of a massive two-loop three-point function. In addition, the program can produce numerical 
results for more general parameter integrals, as they occur for example in phase space integrals 
for multi-particle production with several unresolved massless particles, and offers the possibility 
to include symbolic functions which can be used to define measurement functions like e.g. jet 
algorithms at a later stage. We are looking forward to applications of the program for the 
calculation of higher order corrections to various observables. 
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